GridConvertInteger Subroutine

SUBROUTINE GridConvertInteger &
!
(GridIn, GridOut, cellsize)

USE GeoLib, ONLY: &
!Imported type definitions:
Coordinate, &
! Imported routines:
SetCoord, Convert, ASSIGNMENT(=)


!arguments with intent(in):
TYPE (grid_integer), INTENT (IN) :: GridIn
REAL (KIND = float), OPTIONAL, INTENT (IN) :: cellsize

!arguments with intent (inout):
TYPE (grid_integer), INTENT (INOUT) :: GridOut

!local variables:
TYPE (Coordinate) :: Ain, Bin, Cin, Din !!corner points of input grid
TYPE (Coordinate) :: Aout, Bout, Cout, Dout !!corner points converted to new CRS
TYPE (Coordinate) :: Arec, Brec, Crec, Drec !!corner points of rectangle
TYPE (Coordinate) :: pointNew !!generic point in the new CRS
TYPE (Coordinate) :: pointOld !!generic point in the old CRS
REAL (KIND = float) :: Xdim, Ydim
INTEGER :: newXsize, newYsize
REAL (KIND = float) :: CellSizeX, CellSizeY, newCellSize
INTEGER (KIND = short)  :: ios !!error return code
INTEGER (KIND = short)  :: i, j, iOld, jOld
REAL (KIND = float) :: X, Y
LOGICAL :: checkBound

!------------end of declaration------------------------------------------------

!------------------------------------------------------------------------------
!            calculate corner points
!------------------------------------------------------------------------------
!Set coordinate of corner Ain
CALL SetCoord (GridIn % xllcorner , GridIn % yllcorner + &
               GridIn % idim * GridIn % cellsize, Ain)
!copy corrdinate reference system parameters
Ain % system = GridIn % grid_mapping
Aout %system = GridOut % grid_mapping
!calculate coordinate in the new reference system
CALL Convert (Ain, Aout)

!Set coordinate of corner Bin
CALL SetCoord (GridIn % xllcorner + GridIn % jdim * GridIn % cellsize , &
               GridIn % yllcorner + GridIn % idim * GridIn % cellsize, Bin)
!copy corrdinate reference system parameters
Bin % system = GridIn % grid_mapping
Bout %system = GridOut % grid_mapping
!calculate coordinate in the new reference system
CALL Convert (Bin, Bout)

!Set coordinate of corner Cin
CALL SetCoord (GridIn % xllcorner, GridIn % yllcorner, Cin)
!copy corrdinate reference system parameters
Cin % system = GridIn % grid_mapping
Cout %system = GridOut % grid_mapping
!calculate coordinate in the new reference system
CALL Convert (Cin, Cout)

!Set coordinate of corner Din
CALL SetCoord (GridIn % xllcorner + GridIn % jdim * GridIn % cellsize, &
                GridIn % yllcorner, Din)
!copy corrdinate reference system parameters
Din % system = GridIn % grid_mapping
Dout % system = GridOut % grid_mapping
!calculate coordinate in the new reference system
CALL Convert (Din, Dout)

!------------------------------------------------------------------------------
!        define minimum rectangular view that contains the 4 converted points
!------------------------------------------------------------------------------
!calculate Arec
CALL SetCoord ( MIN (Aout % easting, Cout % easting),  &
                  MAX (Aout % northing, Bout % northing), Arec )
!copy corrdinate reference system parameters
Arec % system = Aout % system

!calculate Brec
CALL SetCoord ( MAX (Bout % easting, Dout % easting),  &
                  MAX (Aout % northing, Bout % northing), Brec )
!copy corrdinate reference system parameters
Brec % system = Bout % system

!calculate Crec
CALL SetCoord ( MIN (Aout % easting, Cout % easting),  &
                  MIN (Cout % northing, Dout % northing), Crec )
!copy corrdinate reference system parameters
Crec % system = Cout % system

!calculate Drec
CALL SetCoord ( MAX (Bout % easting, Dout % easting),  &
                  MIN (Cout % northing, Dout % northing), Drec )
!copy corrdinate reference system parameters
Drec % system = Dout % system


!------------------------------------------------------------------------------
!            calculate cellsize of gridout, number of rows and columns
!------------------------------------------------------------------------------
Xdim = Brec % easting - Arec % easting
Ydim = Brec % northing - Drec % northing
IF (PRESENT (cellsize) ) THEN
  !set cellsize of new grid to cellsize
    newCellSize = cellsize
  !calculate number of columns of new grid
    newXsize = INT ( Xdim / newCellSize ) + 1
  !calculate number of rows of new grid
    newYsize = INT ( Ydim / newCellSize ) + 1
ELSE   
    
    CellSizeX = Xdim / GridIn % jdim

    CellSizeY = Ydim / GridIn % idim

    !set cellsize of new grid to CellSizeX
    newCellSize = CellSizeX 

    !set number of columns of new grid equals to GridIn
    newXsize = GridIn % jdim

    !calculate number of rows of new grid
    newYsize = INT ( Ydim / newCellSize ) + 1
END IF

!------------------------------------------------------------------------------
!                 define new grid
!------------------------------------------------------------------------------
GridOut % jdim = newXsize
GridOut % idim = newYsize
GridOut % standard_name = GridIn % standard_name
GridOut % long_name = GridIn % long_name
GridOut % units = GridIn % units
GridOut % varying_mode = GridIn % varying_mode
GridOut % nodata = GridIn % nodata
GridOut % valid_min = GridIn % valid_min
GridOut % valid_max = GridIn % valid_max
GridOut % reference_time = GridIn % reference_time
GridOut % current_time = GridIn % current_time
GridOut % time_index = GridIn % time_index
GridOut % time_unit = GridIn % time_unit
GridOut % cellsize =  newCellSize
GridOut % xllcorner = Crec % easting
GridOut % yllcorner = Crec % northing
!esri_pe_string !!used by ArcMap 9.2 
IF (ALLOCATED(GridOut % mat)) THEN
  DEALLOCATE (GridOut % mat)
END IF
ALLOCATE ( GridOut % mat ( newYsize, newXsize ), STAT = ios )
IF (ios /= 0) THEN
  CALL Catch ('error', 'GridOperations',  &
  'memory allocation ',  &
  code = memAllocError,argument = 'converted grid' )
ENDIF   

DO i = 1, GridOut % idim
  DO j = 1, GridOut % jdim
    GridOut % mat (i,j) = GridOut % nodata
  END DO
END DO                                 

!------------------------------------------------------------------------------
!                 fill in new grid with values
!------------------------------------------------------------------------------
!set CRS of pointOld and pointNew
pointOld % system = GridIn % grid_mapping
pointNew % system = GridOut % grid_mapping
!loop
DO i = 1, GridOut % idim
  DO j = 1, GridOut % jdim
    !set pointNew
    CALL GetXY (i, j, GridOut, X, Y, checkBound)
    CALL SetCoord (X, Y, pointNew)
    !calculate pointOld
    CALL Convert (pointNew, pointOld)
    CALL GetIJ (pointOld % easting, pointOld % northing, GridIn, iOld, jOld, checkBound)
    IF (checkBound) THEN 
      IF (GridIn % mat (iOld, jOld) /= GridIn % nodata) THEN
        GridOut % mat (i,j) = GridIn % mat (iOld, jOld) 
      END IF    
    END IF
  END DO
END DO     

END SUBROUTINE GridConvertInteger